Robust HMMs
Code to generate state vectors
library(rstan)
options(mc.cores=4)
rstan_options("auto_write" = TRUE)
library(tidyverse)
select_new_state <- function(transition_probs) {
draw <- rmultinom(1, 1, transition_probs)
which(draw==1)
}
simulate_hmm_states <- function(n_steps, state_initial, m_transition_probs) {
states <- vector(length = n_steps)
state <- state_initial
for(i in 1:n_steps) {
transition_probs <- m_transition_probs[state, ]
state <- select_new_state(transition_probs)
states[i] <- state
}
states
}
optimise_repeat <- function(n_opt, model, data_stan) {
best_val <- -Inf
for(i in 1:n_opt) {
fit <- optimizing(model, data=data_stan, as_vector=FALSE)
val <- fit$value
if(val > best_val) {
best_val <- val
fit_best <- fit
}
}
fit_best
}Generate observations from a (1 - w) * normal + w * student_t
r_student_t <- function(mu, sigma, df) {
mu + sigma * rt(n = 1, df = df) #/ sqrt(df / (df - 2))
}
tainted_normal <- function(w, mu, sigma, df) {
u <- runif(1)
if(u < w)
x <- r_student_t(mu, sigma, df)
else
x <- rnorm(1, mu, sigma)
x
}
generate_emissions <- function(states, w, mus, sigmas, dfs) {
xs <- vector(length = length(states))
for(i in seq_along(xs))
xs[i] <- tainted_normal(
w, mus[states[i]], sigmas[states[i]], dfs[states[i]])
xs
}
simulate_hmm <- function(n_steps,
state_initial,
m_transition_probs,
w, mus, sigmas, dfs) {
states <- simulate_hmm_states(n_steps, state_initial, m_transition_probs)
x <- generate_emissions(states, w, mus, sigmas, dfs)
df <- tibble(state=states, obs=x) %>%
mutate(time=seq_along(state)) %>%
mutate(state=as.factor(state))
df
}
mus <- c(1, 8)
sigmas <- c(1, 1)
dfs <- c(3, 3)
K <- 2
m_transition_probs <- matrix(c(0.9, 0.1, 0.1, 0.9), ncol = 2)
n_steps <- 1000
df <- simulate_hmm(n_steps, 1, m_transition_probs, 0, mus, sigmas, dfs)
df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=state))
g <- df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point() +
xlab("Time") +
ylab(latex2exp::TeX("$\\Y_t"))
ggsave("../figures/simulated_normal.pdf", g, width = 6,
height = 3)Fit models
model_normal <- stan_model("hmm.stan")
model_student_t <- stan_model("hmm_student_t.stan")
model_robust <- stan_model("hmm_beta_divergence.stan")
data_stan <- list(
N=nrow(df),
dist=df$obs,
K=2,
beta=0.1
)
n_opt <- 2
fit_normal <- optimise_repeat(n_opt, model_normal, data=data_stan)
fit_student_t <- optimise_repeat(n_opt, model_student_t, data=data_stan)Warning in .local(object, ...) : non-zero return code in optimizingfit_robust <- optimise_repeat(n_opt, model_robust, data=data_stan)
df <- df %>%
mutate(normal=fit_normal$par$state,
student_t=fit_student_t$par$state,
robust=fit_robust$par$state) %>%
mutate(normal=as.factor(normal),
student_t=as.factor(student_t),
robust=as.factor(robust))
df %>%
pivot_longer(-c("time", "obs")) %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=value)) +
facet_wrap(~name)
# some issue with algo for t=1
df$normal[1] <- 1
g <- df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=normal)) +
xlab("Time") +
ylab(latex2exp::TeX("$\\Y_t")) +
scale_color_brewer("State", palette = "Dark2")
ggsave("../figures/simulated_normal_state.pdf", g, width = 6,
height = 3)Now corrupt data
data_stan <- list(
N=nrow(df),
dist=df$obs,
K=2,
beta=0.1
)
n_opt <- 2
fit_normal <- optimise_repeat(n_opt, model_normal, data=data_stan)
fit_student_t <- optimise_repeat(n_opt, model_student_t, data=data_stan)
fit_robust <- optimise_repeat(n_opt, model_robust, data=data_stan)
df <- df %>%
mutate(normal=fit_normal$par$state,
student_t=fit_student_t$par$state,
robust=fit_robust$par$state) %>%
mutate(normal=as.factor(normal),
student_t=as.factor(student_t),
robust=as.factor(robust))# some issue with algo for t=1
df$normal[1] <- 1
g <- df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=normal)) +
xlab("Time") +
ylab(latex2exp::TeX("$\\Y_t")) +
scale_color_brewer("State", palette = "Dark2")
ggsave("../figures/simulated_corrupt_normal_state.pdf", g, width = 6,
height = 3)With 3 states
data_stan <- list(
N=nrow(df),
dist=df$obs,
K=3,
beta=0.1
)
n_opt <- 2
fit_normal <- optimise_repeat(n_opt, model_normal, data=data_stan)
df <- df %>%
mutate(normal=fit_normal$par$state,
student_t=fit_student_t$par$state,
robust=fit_robust$par$state) %>%
mutate(normal=as.factor(normal),
student_t=as.factor(student_t),
robust=as.factor(robust))
# some issue with algo for t=1
df$normal[1] <- 1
g <- df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=normal)) +
xlab("Time") +
ylab(latex2exp::TeX("$\\Y_t")) +
scale_color_brewer("State", palette = "Dark2")
g
ggsave("../figures/simulated_corrupt_normal_state_3.pdf", g, width = 6,
height = 3)
How many states?
find_optimal_k <- function(n_opt, base_parameters, model, data_stan, ks=2:4) {
bics <- vector(length = length(ks))
for(i in seq_along(bics)){
K <- ks[i]
data_stan$K <- K
fit <- optimise_repeat(n_opt, model, data=data_stan)
log_like <- fit$par$log_p;
num_parameters <- K * base_parameters + K * (K - 1)
bics[i] <- 2 * log_like - num_parameters * log(data_stan$N)
}
bics
}
ks <- seq(2, 4, 1)
bics <- find_optimal_k(5, 3, model_normal, data_stan, ks=ks)
plot(ks, bics)
Look at purity of state estimation
mean(df$state==df$normal)[1] 0.0409mean(df$state==df$robust)[1] 0.9991mean(df$state==df$student_t)[1] 0.999Try fitting a DP to fake HMM data
dfs <- c(3, 3)
df <- simulate_hmm(n_steps, 1, m_transition_probs, 0.2, mus, sigmas, dfs)
df %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=state))
dp <- DirichletProcessGaussian(scale(df$obs), g0Priors = c(0, 1, 2, 1),
alphaPriors = c(0.1, 4))
n_iter <- 400
dp_1 <- Fit(dp, n_iter)
|
| | 0%
|
| | 1%
|
|- | 1%
|
|- | 2%
|
|- | 3%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 5%
|
|--- | 6%
|
|--- | 7%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 9%
|
|----- | 10%
|
|----- | 11%
|
|------ | 11%
|
|------ | 12%
|
|------ | 13%
|
|------- | 13%
|
|------- | 14%
|
|------- | 15%
|
|-------- | 15%
|
|-------- | 16%
|
|-------- | 17%
|
|--------- | 17%
|
|--------- | 18%
|
|--------- | 19%
|
|---------- | 19%
|
|---------- | 20%
|
|---------- | 21%
|
|----------- | 21%
|
|----------- | 22%
|
|----------- | 23%
|
|------------ | 23%
|
|------------ | 24%
|
|------------ | 25%
|
|------------- | 25%
|
|------------- | 26%
|
|------------- | 27%
|
|-------------- | 27%
|
|-------------- | 28%
|
|-------------- | 29%
|
|--------------- | 29%
|
|--------------- | 30%
|
|--------------- | 31%
|
|---------------- | 31%
|
|---------------- | 32%
|
|---------------- | 33%
|
|----------------- | 33%
|
|----------------- | 34%
|
|----------------- | 35%
|
|------------------ | 35%
|
|------------------ | 36%
|
|------------------ | 37%
|
|------------------- | 37%
|
|------------------- | 38%
|
|------------------- | 39%
|
|-------------------- | 39%
|
|-------------------- | 40%
|
|-------------------- | 41%
|
|--------------------- | 41%
|
|--------------------- | 42%
|
|--------------------- | 43%
|
|---------------------- | 43%
|
|---------------------- | 44%
|
|---------------------- | 45%
|
|----------------------- | 45%
|
|----------------------- | 46%
|
|----------------------- | 47%
|
|------------------------ | 47%
|
|------------------------ | 48%
|
|------------------------ | 49%
|
|------------------------- | 49%
|
|------------------------- | 50%
|
|------------------------- | 51%
|
|-------------------------- | 51%
|
|-------------------------- | 52%
|
|-------------------------- | 53%
|
|--------------------------- | 53%
|
|--------------------------- | 54%
|
|--------------------------- | 55%
|
|---------------------------- | 55%
|
|---------------------------- | 56%
|
|---------------------------- | 57%
|
|----------------------------- | 57%
|
|----------------------------- | 58%
|
|----------------------------- | 59%
|
|------------------------------ | 59%
|
|------------------------------ | 60%
|
|------------------------------ | 61%
|
|------------------------------- | 61%
|
|------------------------------- | 62%
|
|------------------------------- | 63%
|
|-------------------------------- | 63%
|
|-------------------------------- | 64%
|
|-------------------------------- | 65%
|
|--------------------------------- | 65%
|
|--------------------------------- | 66%
|
|--------------------------------- | 67%
|
|---------------------------------- | 67%
|
|---------------------------------- | 68%
|
|---------------------------------- | 69%
|
|----------------------------------- | 69%
|
|----------------------------------- | 70%
|
|----------------------------------- | 71%
|
|------------------------------------ | 71%
|
|------------------------------------ | 72%
|
|------------------------------------ | 73%
|
|------------------------------------- | 73%
|
|------------------------------------- | 74%
|
|------------------------------------- | 75%
|
|-------------------------------------- | 75%
|
|-------------------------------------- | 76%
|
|-------------------------------------- | 77%
|
|--------------------------------------- | 77%
|
|--------------------------------------- | 78%
|
|--------------------------------------- | 79%
|
|---------------------------------------- | 79%
|
|---------------------------------------- | 80%
|
|---------------------------------------- | 81%
|
|----------------------------------------- | 81%
|
|----------------------------------------- | 82%
|
|----------------------------------------- | 83%
|
|------------------------------------------ | 83%
|
|------------------------------------------ | 84%
|
|------------------------------------------ | 85%
|
|------------------------------------------- | 85%
|
|------------------------------------------- | 86%
|
|------------------------------------------- | 87%
|
|-------------------------------------------- | 87%
|
|-------------------------------------------- | 88%
|
|-------------------------------------------- | 89%
|
|--------------------------------------------- | 89%
|
|--------------------------------------------- | 90%
|
|--------------------------------------------- | 91%
|
|---------------------------------------------- | 91%
|
|---------------------------------------------- | 92%
|
|---------------------------------------------- | 93%
|
|----------------------------------------------- | 93%
|
|----------------------------------------------- | 94%
|
|----------------------------------------------- | 95%
|
|------------------------------------------------ | 95%
|
|------------------------------------------------ | 96%
|
|------------------------------------------------ | 97%
|
|------------------------------------------------- | 97%
|
|------------------------------------------------- | 98%
|
|------------------------------------------------- | 99%
|
|--------------------------------------------------| 99%
|
|--------------------------------------------------| 100%dp_2 <- Fit(dp, n_iter)
|
| | 0%
|
| | 1%
|
|- | 1%
|
|- | 2%
|
|- | 3%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 5%
|
|--- | 6%
|
|--- | 7%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 9%
|
|----- | 10%
|
|----- | 11%
|
|------ | 11%
|
|------ | 12%
|
|------ | 13%
|
|------- | 13%
|
|------- | 14%
|
|------- | 15%
|
|-------- | 15%
|
|-------- | 16%
|
|-------- | 17%
|
|--------- | 17%
|
|--------- | 18%
|
|--------- | 19%
|
|---------- | 19%
|
|---------- | 20%
|
|---------- | 21%
|
|----------- | 21%
|
|----------- | 22%
|
|----------- | 23%
|
|------------ | 23%
|
|------------ | 24%
|
|------------ | 25%
|
|------------- | 25%
|
|------------- | 26%
|
|------------- | 27%
|
|-------------- | 27%
|
|-------------- | 28%
|
|-------------- | 29%
|
|--------------- | 29%
|
|--------------- | 30%
|
|--------------- | 31%
|
|---------------- | 31%
|
|---------------- | 32%
|
|---------------- | 33%
|
|----------------- | 33%
|
|----------------- | 34%
|
|----------------- | 35%
|
|------------------ | 35%
|
|------------------ | 36%
|
|------------------ | 37%
|
|------------------- | 37%
|
|------------------- | 38%
|
|------------------- | 39%
|
|-------------------- | 39%
|
|-------------------- | 40%
|
|-------------------- | 41%
|
|--------------------- | 41%
|
|--------------------- | 42%
|
|--------------------- | 43%
|
|---------------------- | 43%
|
|---------------------- | 44%
|
|---------------------- | 45%
|
|----------------------- | 45%
|
|----------------------- | 46%
|
|----------------------- | 47%
|
|------------------------ | 47%
|
|------------------------ | 48%
|
|------------------------ | 49%
|
|------------------------- | 49%
|
|------------------------- | 50%
|
|------------------------- | 51%
|
|-------------------------- | 51%
|
|-------------------------- | 52%
|
|-------------------------- | 53%
|
|--------------------------- | 53%
|
|--------------------------- | 54%
|
|--------------------------- | 55%
|
|---------------------------- | 55%
|
|---------------------------- | 56%
|
|---------------------------- | 57%
|
|----------------------------- | 57%
|
|----------------------------- | 58%
|
|----------------------------- | 59%
|
|------------------------------ | 59%
|
|------------------------------ | 60%
|
|------------------------------ | 61%
|
|------------------------------- | 61%
|
|------------------------------- | 62%
|
|------------------------------- | 63%
|
|-------------------------------- | 63%
|
|-------------------------------- | 64%
|
|-------------------------------- | 65%
|
|--------------------------------- | 65%
|
|--------------------------------- | 66%
|
|--------------------------------- | 67%
|
|---------------------------------- | 67%
|
|---------------------------------- | 68%
|
|---------------------------------- | 69%
|
|----------------------------------- | 69%
|
|----------------------------------- | 70%
|
|----------------------------------- | 71%
|
|------------------------------------ | 71%
|
|------------------------------------ | 72%
|
|------------------------------------ | 73%
|
|------------------------------------- | 73%
|
|------------------------------------- | 74%
|
|------------------------------------- | 75%
|
|-------------------------------------- | 75%
|
|-------------------------------------- | 76%
|
|-------------------------------------- | 77%
|
|--------------------------------------- | 77%
|
|--------------------------------------- | 78%
|
|--------------------------------------- | 79%
|
|---------------------------------------- | 79%
|
|---------------------------------------- | 80%
|
|---------------------------------------- | 81%
|
|----------------------------------------- | 81%
|
|----------------------------------------- | 82%
|
|----------------------------------------- | 83%
|
|------------------------------------------ | 83%
|
|------------------------------------------ | 84%
|
|------------------------------------------ | 85%
|
|------------------------------------------- | 85%
|
|------------------------------------------- | 86%
|
|------------------------------------------- | 87%
|
|-------------------------------------------- | 87%
|
|-------------------------------------------- | 88%
|
|-------------------------------------------- | 89%
|
|--------------------------------------------- | 89%
|
|--------------------------------------------- | 90%
|
|--------------------------------------------- | 91%
|
|---------------------------------------------- | 91%
|
|---------------------------------------------- | 92%
|
|---------------------------------------------- | 93%
|
|----------------------------------------------- | 93%
|
|----------------------------------------------- | 94%
|
|----------------------------------------------- | 95%
|
|------------------------------------------------ | 95%
|
|------------------------------------------------ | 96%
|
|------------------------------------------------ | 97%
|
|------------------------------------------------- | 97%
|
|------------------------------------------------- | 98%
|
|------------------------------------------------- | 99%
|
|--------------------------------------------------| 99%
|
|--------------------------------------------------| 100%n_clusters_1 <- map_dbl(dp_1$labelsChain, ~max(.))
n_clusters_2 <- map_dbl(dp_2$labelsChain, ~max(.))
tibble(n1=n_clusters_1, n2=n_clusters_2) %>%
mutate(iter=seq_along(n1)) %>%
pivot_longer(-iter) %>%
ggplot(aes(x=iter, y=value, colour=name)) +
geom_line()
dp_1$numberClusters[1] 7dp_2$numberClusters[1] 1dp_1$alpha[1] 0.455211df %>%
mutate(cluster=dp_2$clusterLabels) %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=cluster))
Try tempering the likelihood
Likelihood.tempered <- function(mdobj, x, theta){
beta <- 1
val1 <- 1 / beta * dnorm(x, theta[[1]], theta[[2]])^beta
# val2 <- (1 / (beta + 1)) * integrate(function(x) dnorm(x, theta[[1]], theta[[2]])^(beta + 1), -20, 20)$value
as.numeric(val1)
}
PriorDraw.tempered <- function(mdobj, n=1){
theta <- list()
theta[[1]] = array(rnorm(n, mdobj$priorParameters[1], mdobj$priorParameters[2]), dim=c(1, 1, n))
theta[[2]] = array(rexp(n, mdobj$priorParameters[3]), dim=c(1, 1, n))
theta
}
PriorDensity.tempered <- function(mdobj, theta){
priorParameters <- mdobj$priorParameters
thetaDensity <- dnorm(theta[[1]], priorParameters[1], priorParameters[2])
thetaDensity <- thetaDensity * dexp(theta[[2]], priorParameters[3])
as.numeric(thetaDensity)
}
MhParameterProposal.tempered <- function(mdobj, oldParams){
mhStepSize <- mdobj$mhStepSize
newParams <- oldParams
newParams[[1]] <- oldParams[[1]] + mhStepSize[1]*rnorm(1)
newParams[[2]] <- abs(oldParams[[2]] + mhStepSize[2]*rnorm(1))
newParams
}
temperedMd <- MixingDistribution(distribution = "tempered",
priorParameters = c(0, 5, 0.01),
conjugate = "nonconjugate",
mhStepSize = c(1, 0.1)
)
dp <- DirichletProcessCreate(scale(df$obs), temperedMd)
dp <- Initialise(dp, numInitialClusters = 3)Accept Ratio: 0.695 dp <- Fit(dp, 200)
|
| | 0%
|
| | 1%
|
|- | 2%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 6%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 10%
|
|------ | 11%
|
|------ | 12%
|
|------ | 13%
|
|------- | 14%
|
|-------- | 15%
|
|-------- | 16%
|
|-------- | 17%
|
|--------- | 18%
|
|---------- | 19%
|
|---------- | 20%
|
|---------- | 21%
|
|----------- | 22%
|
|------------ | 23%
|
|------------ | 24%
|
|------------ | 25%
|
|------------- | 26%
|
|-------------- | 27%
|
|-------------- | 28%
|
|-------------- | 29%
|
|--------------- | 30%
|
|---------------- | 31%
|
|---------------- | 32%
|
|---------------- | 33%
|
|----------------- | 34%
|
|------------------ | 35%
|
|------------------ | 36%
|
|------------------ | 37%
|
|------------------- | 38%
|
|-------------------- | 39%
|
|-------------------- | 40%
|
|-------------------- | 41%
|
|--------------------- | 42%
|
|---------------------- | 43%
|
|---------------------- | 44%
|
|---------------------- | 45%
|
|----------------------- | 46%
|
|------------------------ | 47%
|
|------------------------ | 48%
|
|------------------------ | 49%
|
|------------------------- | 50%
|
|-------------------------- | 51%
|
|-------------------------- | 52%
|
|-------------------------- | 53%
|
|--------------------------- | 54%
|
|---------------------------- | 55%
|
|---------------------------- | 56%
|
|---------------------------- | 57%
|
|----------------------------- | 58%
|
|------------------------------ | 59%
|
|------------------------------ | 60%
|
|------------------------------ | 61%
|
|------------------------------- | 62%
|
|-------------------------------- | 63%
|
|-------------------------------- | 64%
|
|-------------------------------- | 65%
|
|--------------------------------- | 66%
|
|---------------------------------- | 67%
|
|---------------------------------- | 68%
|
|---------------------------------- | 69%
|
|----------------------------------- | 70%
|
|------------------------------------ | 71%
|
|------------------------------------ | 72%
|
|------------------------------------ | 73%
|
|------------------------------------- | 74%
|
|-------------------------------------- | 75%
|
|-------------------------------------- | 76%
|
|-------------------------------------- | 77%
|
|--------------------------------------- | 78%
|
|---------------------------------------- | 79%
|
|---------------------------------------- | 80%
|
|---------------------------------------- | 81%
|
|----------------------------------------- | 82%
|
|------------------------------------------ | 83%
|
|------------------------------------------ | 84%
|
|------------------------------------------ | 85%
|
|------------------------------------------- | 86%
|
|-------------------------------------------- | 87%
|
|-------------------------------------------- | 88%
|
|-------------------------------------------- | 89%
|
|--------------------------------------------- | 90%
|
|---------------------------------------------- | 91%
|
|---------------------------------------------- | 92%
|
|---------------------------------------------- | 93%
|
|----------------------------------------------- | 94%
|
|------------------------------------------------ | 95%
|
|------------------------------------------------ | 96%
|
|------------------------------------------------ | 97%
|
|------------------------------------------------- | 98%
|
|--------------------------------------------------| 99%
|
|--------------------------------------------------| 100%df %>%
mutate(cluster=dp$clusterLabels) %>%
ggplot(aes(x=time, y=obs)) +
geom_line() +
geom_point(aes(colour=cluster))
Try fitting DP to fake 2 mean data
Real data
Fit DP to data
dp_1$alpha
[1] 0.346678Try with more angle data too
df_s <- df %>%
filter(id==1) %>%
slice(1:200) %>%
select(Dist, Angle) %>%
scale()
df_s %>%
as.data.frame() %>%
mutate(time=seq_along(Dist)) %>%
pivot_longer(-time) %>%
ggplot(aes(x=time, y=value)) +
geom_line() +
facet_wrap(~name)
dp_1 <- DirichletProcessMvnormal(df_s, numInitialClusters = 2, alphaPriors = c(0.001, 4))
dp_2 <- DirichletProcessMvnormal(df_s, numInitialClusters = 2, alphaPriors = c(0.001, 4))
n_iter <- 200
dp_1 <- Fit(dp_1, n_iter)
|
| | 0%
|
| | 1%
|
|- | 2%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 6%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 10%
|
|------ | 11%
|
|------ | 12%
|
|------ | 13%
|
|------- | 14%
|
|-------- | 15%
|
|-------- | 16%
|
|-------- | 17%
|
|--------- | 18%
|
|---------- | 19%
|
|---------- | 20%
|
|---------- | 21%
|
|----------- | 22%
|
|------------ | 23%
|
|------------ | 24%
|
|------------ | 25%
|
|------------- | 26%
|
|-------------- | 27%
|
|-------------- | 28%
|
|-------------- | 29%
|
|--------------- | 30%
|
|---------------- | 31%
|
|---------------- | 32%
|
|---------------- | 33%
|
|----------------- | 34%
|
|------------------ | 35%
|
|------------------ | 36%
|
|------------------ | 37%
|
|------------------- | 38%
|
|-------------------- | 39%
|
|-------------------- | 40%
|
|-------------------- | 41%
|
|--------------------- | 42%
|
|---------------------- | 43%
|
|---------------------- | 44%
|
|---------------------- | 45%
|
|----------------------- | 46%
|
|------------------------ | 47%
|
|------------------------ | 48%
|
|------------------------ | 49%
|
|------------------------- | 50%
|
|-------------------------- | 51%
|
|-------------------------- | 52%
|
|-------------------------- | 53%
|
|--------------------------- | 54%
|
|---------------------------- | 55%
|
|---------------------------- | 56%
|
|---------------------------- | 57%
|
|----------------------------- | 58%
|
|------------------------------ | 59%
|
|------------------------------ | 60%
|
|------------------------------ | 61%
|
|------------------------------- | 62%
|
|-------------------------------- | 63%
|
|-------------------------------- | 64%
|
|-------------------------------- | 65%
|
|--------------------------------- | 66%
|
|---------------------------------- | 67%
|
|---------------------------------- | 68%
|
|---------------------------------- | 69%
|
|----------------------------------- | 70%
|
|------------------------------------ | 71%
|
|------------------------------------ | 72%
|
|------------------------------------ | 73%
|
|------------------------------------- | 74%
|
|-------------------------------------- | 75%
|
|-------------------------------------- | 76%
|
|-------------------------------------- | 77%
|
|--------------------------------------- | 78%
|
|---------------------------------------- | 79%
|
|---------------------------------------- | 80%
|
|---------------------------------------- | 81%
|
|----------------------------------------- | 82%
|
|------------------------------------------ | 83%
|
|------------------------------------------ | 84%
|
|------------------------------------------ | 85%
|
|------------------------------------------- | 86%
|
|-------------------------------------------- | 87%
|
|-------------------------------------------- | 88%
|
|-------------------------------------------- | 89%
|
|--------------------------------------------- | 90%
|
|---------------------------------------------- | 91%
|
|---------------------------------------------- | 92%
|
|---------------------------------------------- | 93%
|
|----------------------------------------------- | 94%
|
|------------------------------------------------ | 95%
|
|------------------------------------------------ | 96%
|
|------------------------------------------------ | 97%
|
|------------------------------------------------- | 98%
|
|--------------------------------------------------| 99%
|
|--------------------------------------------------| 100%dp_2 <- Fit(dp_2, n_iter)
|
| | 0%
|
| | 1%
|
|- | 2%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 6%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 10%
|
|------ | 11%
|
|------ | 12%
|
|------ | 13%
|
|------- | 14%
|
|-------- | 15%
|
|-------- | 16%
|
|-------- | 17%
|
|--------- | 18%
|
|---------- | 19%
|
|---------- | 20%
|
|---------- | 21%
|
|----------- | 22%
|
|------------ | 23%
|
|------------ | 24%
|
|------------ | 25%
|
|------------- | 26%
|
|-------------- | 27%
|
|-------------- | 28%
|
|-------------- | 29%
|
|--------------- | 30%
|
|---------------- | 31%
|
|---------------- | 32%
|
|---------------- | 33%
|
|----------------- | 34%
|
|------------------ | 35%
|
|------------------ | 36%
|
|------------------ | 37%
|
|------------------- | 38%
|
|-------------------- | 39%
|
|-------------------- | 40%
|
|-------------------- | 41%
|
|--------------------- | 42%
|
|---------------------- | 43%
|
|---------------------- | 44%
|
|---------------------- | 45%
|
|----------------------- | 46%
|
|------------------------ | 47%
|
|------------------------ | 48%
|
|------------------------ | 49%
|
|------------------------- | 50%
|
|-------------------------- | 51%
|
|-------------------------- | 52%
|
|-------------------------- | 53%
|
|--------------------------- | 54%
|
|---------------------------- | 55%
|
|---------------------------- | 56%
|
|---------------------------- | 57%
|
|----------------------------- | 58%
|
|------------------------------ | 59%
|
|------------------------------ | 60%
|
|------------------------------ | 61%
|
|------------------------------- | 62%
|
|-------------------------------- | 63%
|
|-------------------------------- | 64%
|
|-------------------------------- | 65%
|
|--------------------------------- | 66%
|
|---------------------------------- | 67%
|
|---------------------------------- | 68%
|
|---------------------------------- | 69%
|
|----------------------------------- | 70%
|
|------------------------------------ | 71%
|
|------------------------------------ | 72%
|
|------------------------------------ | 73%
|
|------------------------------------- | 74%
|
|-------------------------------------- | 75%
|
|-------------------------------------- | 76%
|
|-------------------------------------- | 77%
|
|--------------------------------------- | 78%
|
|---------------------------------------- | 79%
|
|---------------------------------------- | 80%
|
|---------------------------------------- | 81%
|
|----------------------------------------- | 82%
|
|------------------------------------------ | 83%
|
|------------------------------------------ | 84%
|
|------------------------------------------ | 85%
|
|------------------------------------------- | 86%
|
|-------------------------------------------- | 87%
|
|-------------------------------------------- | 88%
|
|-------------------------------------------- | 89%
|
|--------------------------------------------- | 90%
|
|---------------------------------------------- | 91%
|
|---------------------------------------------- | 92%
|
|---------------------------------------------- | 93%
|
|----------------------------------------------- | 94%
|
|------------------------------------------------ | 95%
|
|------------------------------------------------ | 96%
|
|------------------------------------------------ | 97%
|
|------------------------------------------------- | 98%
|
|--------------------------------------------------| 99%
|
|--------------------------------------------------| 100%n_clusters_1 <- map_dbl(dp_1$labelsChain, ~max(.))
n_clusters_2 <- map_dbl(dp_2$labelsChain, ~max(.))
tibble(n1=n_clusters_1, n2=n_clusters_2) %>%
mutate(iter=seq_along(n1)) %>%
pivot_longer(-iter) %>%
ggplot(aes(x=iter, y=value, colour=name)) +
geom_line()
Look at clusters
DPM on Mathematica HMM data
x <- read.csv("../data/hmm_test.csv", header = FALSE)
x <- x %>%
slice(1:nrow(x))
dp_1 <- DirichletProcessGaussian(
scale(x$V2), alphaPriors = c(4, 4))
dp_2 <- DirichletProcessGaussian(
scale(x$V2), alphaPriors = c(4, 4))
dp_1 <- Initialise(dp_1, numInitialClusters = 1)
dp_2 <- Initialise(dp_2, numInitialClusters = 10)
n_iter <- 2000
dp_1 <- Fit(dp, n_iter)
|
| | 0%
|
| | 1%
|
|- | 1%
|
|- | 2%
|
|- | 3%
|
|-- | 3%
|
|-- | 4%
|
|-- | 5%
|
|--- | 5%
|
|--- | 6%
|
|--- | 7%
|
|---- | 7%
|
|---- | 8%
|
|---- | 9%
|
|----- | 9%
|
|----- | 10%
|
|----- | 11%
|
|------ | 11%
|
|------ | 12%LS0tCnRpdGxlOiAiUm9idXN0IEhNTXMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCkNvZGUgdG8gZ2VuZXJhdGUgc3RhdGUgdmVjdG9ycwpgYGB7cn0Kc291cmNlKCJoZWxwZXIuUiIpCmBgYAoKR2VuZXJhdGUgb2JzZXJ2YXRpb25zIGZyb20gYSAoMSAtIHcpICogbm9ybWFsICsgdyAqIHN0dWRlbnRfdApgYGB7cn0KbXVzIDwtIGMoMSwgOCkKc2lnbWFzIDwtIGMoMSwgMSkKZGZzIDwtIGMoMywgMykKSyA8LSAyCm1fdHJhbnNpdGlvbl9wcm9icyA8LSBtYXRyaXgoYygwLjksIDAuMSwgMC4xLCAwLjkpLCBuY29sID0gMikKCm5fc3RlcHMgPC0gMTAwMApkZiA8LSBzaW11bGF0ZV9obW0obl9zdGVwcywgMSwgbV90cmFuc2l0aW9uX3Byb2JzLCAwLCBtdXMsIHNpZ21hcywgZGZzKQoKZGYgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPXN0YXRlKSkKYGBgCgpgYGB7cn0KZyA8LSBkZiAlPiUgCiAgZ2dwbG90KGFlcyh4PXRpbWUsIHk9b2JzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KCkgKwogIHhsYWIoIlRpbWUiKSArCiAgeWxhYihsYXRleDJleHA6OlRlWCgiJFxcWV90IikpCgpnZ3NhdmUoIi4uL2ZpZ3VyZXMvc2ltdWxhdGVkX25vcm1hbC5wZGYiLCBnLCB3aWR0aCA9IDYsCiAgICAgICBoZWlnaHQgPSAzKQpgYGAKCgpGaXQgbW9kZWxzCmBgYHtyfQptb2RlbF9ub3JtYWwgPC0gc3Rhbl9tb2RlbCgiaG1tLnN0YW4iKQptb2RlbF9zdHVkZW50X3QgPC0gc3Rhbl9tb2RlbCgiaG1tX3N0dWRlbnRfdC5zdGFuIikKbW9kZWxfcm9idXN0IDwtIHN0YW5fbW9kZWwoImhtbV9iZXRhX2RpdmVyZ2VuY2Uuc3RhbiIpCgpkYXRhX3N0YW4gPC0gbGlzdCgKICBOPW5yb3coZGYpLAogIGRpc3Q9ZGYkb2JzLAogIEs9MiwKICBiZXRhPTAuMQogICkKCm5fb3B0IDwtIDIKZml0X25vcm1hbCA8LSBvcHRpbWlzZV9yZXBlYXQobl9vcHQsIG1vZGVsX25vcm1hbCwgZGF0YT1kYXRhX3N0YW4pCmZpdF9zdHVkZW50X3QgPC0gb3B0aW1pc2VfcmVwZWF0KG5fb3B0LCBtb2RlbF9zdHVkZW50X3QsIGRhdGE9ZGF0YV9zdGFuKQpmaXRfcm9idXN0IDwtIG9wdGltaXNlX3JlcGVhdChuX29wdCwgbW9kZWxfcm9idXN0LCBkYXRhPWRhdGFfc3RhbikKCmRmIDwtIGRmICU+JSAKICBtdXRhdGUobm9ybWFsPWZpdF9ub3JtYWwkcGFyJHN0YXRlLAogICAgICAgICBzdHVkZW50X3Q9Zml0X3N0dWRlbnRfdCRwYXIkc3RhdGUsCiAgICAgICAgIHJvYnVzdD1maXRfcm9idXN0JHBhciRzdGF0ZSkgJT4lIAogIG11dGF0ZShub3JtYWw9YXMuZmFjdG9yKG5vcm1hbCksCiAgICAgICAgIHN0dWRlbnRfdD1hcy5mYWN0b3Ioc3R1ZGVudF90KSwKICAgICAgICAgcm9idXN0PWFzLmZhY3Rvcihyb2J1c3QpKQoKZGYgJT4lIAogIHBpdm90X2xvbmdlcigtYygidGltZSIsICJvYnMiKSkgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPXZhbHVlKSkgKwogIGZhY2V0X3dyYXAofm5hbWUpCmBgYAoKYGBge3J9CiMgc29tZSBpc3N1ZSB3aXRoIGFsZ28gZm9yIHQ9MQpkZiRub3JtYWxbMV0gPC0gMQoKZyA8LSBkZiAlPiUKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT1vYnMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1ub3JtYWwpKSArCiAgeGxhYigiVGltZSIpICsKICB5bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3QiKSkgKwogIHNjYWxlX2NvbG9yX2JyZXdlcigiU3RhdGUiLCBwYWxldHRlID0gIkRhcmsyIikKCmdnc2F2ZSgiLi4vZmlndXJlcy9zaW11bGF0ZWRfbm9ybWFsX3N0YXRlLnBkZiIsIGcsIHdpZHRoID0gNiwKICAgICAgIGhlaWdodCA9IDMpCmBgYAoKTm93IGNvcnJ1cHQgZGF0YQpgYGB7cn0KZGZzIDwtIGMoMS41LCAxLjUpCmRmIDwtIHNpbXVsYXRlX2htbShuX3N0ZXBzLCAxLCBtX3RyYW5zaXRpb25fcHJvYnMsIDAuMiwgbXVzLCBzaWdtYXMsIGRmcykKCmcgPC0gZGYgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPXN0YXRlKSkgKwogIHhsYWIoIlRpbWUiKSArCiAgeWxhYihsYXRleDJleHA6OlRlWCgiJFxcWV90IikpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpCmcKYGBgCgpgYGB7cn0KZGF0YV9zdGFuIDwtIGxpc3QoCiAgTj1ucm93KGRmKSwKICBkaXN0PWRmJG9icywKICBLPTIsCiAgYmV0YT0wLjEKICApCgpuX29wdCA8LSAyCmZpdF9ub3JtYWwgPC0gb3B0aW1pc2VfcmVwZWF0KG5fb3B0LCBtb2RlbF9ub3JtYWwsIGRhdGE9ZGF0YV9zdGFuKQpmaXRfc3R1ZGVudF90IDwtIG9wdGltaXNlX3JlcGVhdChuX29wdCwgbW9kZWxfc3R1ZGVudF90LCBkYXRhPWRhdGFfc3RhbikKZml0X3JvYnVzdCA8LSBvcHRpbWlzZV9yZXBlYXQobl9vcHQsIG1vZGVsX3JvYnVzdCwgZGF0YT1kYXRhX3N0YW4pCgpkZiA8LSBkZiAlPiUgCiAgbXV0YXRlKG5vcm1hbD1maXRfbm9ybWFsJHBhciRzdGF0ZSwKICAgICAgICAgc3R1ZGVudF90PWZpdF9zdHVkZW50X3QkcGFyJHN0YXRlLAogICAgICAgICByb2J1c3Q9Zml0X3JvYnVzdCRwYXIkc3RhdGUpICU+JSAKICBtdXRhdGUobm9ybWFsPWFzLmZhY3Rvcihub3JtYWwpLAogICAgICAgICBzdHVkZW50X3Q9YXMuZmFjdG9yKHN0dWRlbnRfdCksCiAgICAgICAgIHJvYnVzdD1hcy5mYWN0b3Iocm9idXN0KSkKYGBgCgpgYGB7cn0KIyBzb21lIGlzc3VlIHdpdGggYWxnbyBmb3IgdD0xCmRmJG5vcm1hbFsxXSA8LSAxCgpnIDwtIGRmICU+JQogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPW5vcm1hbCkpICsKICB4bGFiKCJUaW1lIikgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfdCIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKQpgYGAKCldpdGggMyBzdGF0ZXMKYGBge3J9CmRhdGFfc3RhbiA8LSBsaXN0KAogIE49bnJvdyhkZiksCiAgZGlzdD1kZiRvYnMsCiAgSz0zLAogIGJldGE9MC4xCiAgKQoKbl9vcHQgPC0gMgpmaXRfbm9ybWFsIDwtIG9wdGltaXNlX3JlcGVhdChuX29wdCwgbW9kZWxfbm9ybWFsLCBkYXRhPWRhdGFfc3RhbikKCmRmIDwtIGRmICU+JSAKICBtdXRhdGUobm9ybWFsPWZpdF9ub3JtYWwkcGFyJHN0YXRlLAogICAgICAgICBzdHVkZW50X3Q9Zml0X3N0dWRlbnRfdCRwYXIkc3RhdGUsCiAgICAgICAgIHJvYnVzdD1maXRfcm9idXN0JHBhciRzdGF0ZSkgJT4lIAogIG11dGF0ZShub3JtYWw9YXMuZmFjdG9yKG5vcm1hbCksCiAgICAgICAgIHN0dWRlbnRfdD1hcy5mYWN0b3Ioc3R1ZGVudF90KSwKICAgICAgICAgcm9idXN0PWFzLmZhY3Rvcihyb2J1c3QpKQoKIyBzb21lIGlzc3VlIHdpdGggYWxnbyBmb3IgdD0xCmRmJG5vcm1hbFsxXSA8LSAxCgpnIDwtIGRmICU+JQogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPW5vcm1hbCkpICsKICB4bGFiKCJUaW1lIikgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfdCIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKQpnCmBgYAoKCgpIb3cgbWFueSBzdGF0ZXM/CmBgYHtyfQpmaW5kX29wdGltYWxfayA8LSBmdW5jdGlvbihuX29wdCwgYmFzZV9wYXJhbWV0ZXJzLCBtb2RlbCwgZGF0YV9zdGFuLCBrcz0yOjQpIHsKICBiaWNzIDwtIHZlY3RvcihsZW5ndGggPSBsZW5ndGgoa3MpKQogIGZvcihpIGluIHNlcV9hbG9uZyhiaWNzKSl7CiAgICBLIDwtIGtzW2ldCiAgICBkYXRhX3N0YW4kSyA8LSBLIAogICAgZml0IDwtIG9wdGltaXNlX3JlcGVhdChuX29wdCwgbW9kZWwsIGRhdGE9ZGF0YV9zdGFuKQogICAgbG9nX2xpa2UgPC0gZml0JHBhciRsb2dfcDsKICAgIG51bV9wYXJhbWV0ZXJzIDwtIEsgKiBiYXNlX3BhcmFtZXRlcnMgKyBLICogKEsgLSAxKQogICAgYmljc1tpXSA8LSAyICogbG9nX2xpa2UgLSBudW1fcGFyYW1ldGVycyAqIGxvZyhkYXRhX3N0YW4kTikKICB9CiAgYmljcwp9CgprcyA8LSBzZXEoMiwgNCwgMSkKYmljcyA8LSBmaW5kX29wdGltYWxfayg1LCAzLCBtb2RlbF9ub3JtYWwsIGRhdGFfc3Rhbiwga3M9a3MpCnBsb3Qoa3MsIGJpY3MpCmBgYAoKCkxvb2sgYXQgcHVyaXR5IG9mIHN0YXRlIGVzdGltYXRpb24KYGBge3J9Cm1lYW4oZGYkc3RhdGU9PWRmJG5vcm1hbCkKbWVhbihkZiRzdGF0ZT09ZGYkcm9idXN0KQptZWFuKGRmJHN0YXRlPT1kZiRzdHVkZW50X3QpCmBgYAoKVHJ5IGZpdHRpbmcgYSBEUCB0byBmYWtlIEhNTSBkYXRhCmBgYHtyfQpuX3N0ZXBzIDwtIDIwMDAKZGYgPC0gc2ltdWxhdGVfaG1tKG5fc3RlcHMsIDEsIG1fdHJhbnNpdGlvbl9wcm9icywgMCwgbXVzLCBzaWdtYXMsIGRmcykKCmRmICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT1vYnMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1zdGF0ZSkpCgpsaWJyYXJ5KGRpcmljaGxldHByb2Nlc3MpCmRwIDwtIERpcmljaGxldFByb2Nlc3NHYXVzc2lhbihzY2FsZShkZiRvYnMpLCBnMFByaW9ycyA9IGMoMCwgMSwgMiwgMSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhbHBoYVByaW9ycyA9IGMoMC4xLCA0KSkKbl9pdGVyIDwtIDQwMApkcF8xIDwtIEZpdChkcCwgbl9pdGVyKQpkcF8yIDwtIEZpdChkcCwgbl9pdGVyKQoKbl9jbHVzdGVyc18xIDwtIG1hcF9kYmwoZHBfMSRsYWJlbHNDaGFpbiwgfm1heCguKSkKbl9jbHVzdGVyc18yIDwtIG1hcF9kYmwoZHBfMiRsYWJlbHNDaGFpbiwgfm1heCguKSkKCnRpYmJsZShuMT1uX2NsdXN0ZXJzXzEsIG4yPW5fY2x1c3RlcnNfMikgJT4lIAogIG11dGF0ZShpdGVyPXNlcV9hbG9uZyhuMSkpICU+JSAKICBwaXZvdF9sb25nZXIoLWl0ZXIpICU+JSAKICBnZ3Bsb3QoYWVzKHg9aXRlciwgeT12YWx1ZSwgY29sb3VyPW5hbWUpKSArCiAgZ2VvbV9saW5lKCkKCmRwXzEkbnVtYmVyQ2x1c3RlcnMKZHBfMiRudW1iZXJDbHVzdGVycwpkcF8xJGFscGhhCgpkZiAlPiUgCiAgbXV0YXRlKGNsdXN0ZXI9ZHBfMiRjbHVzdGVyTGFiZWxzKSAlPiUgCiAgZ2dwbG90KGFlcyh4PXRpbWUsIHk9b2JzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvdXI9Y2x1c3RlcikpCmBgYAoKYGBge3J9CmRmcyA8LSBjKDMsIDMpCmRmIDwtIHNpbXVsYXRlX2htbShuX3N0ZXBzLCAxLCBtX3RyYW5zaXRpb25fcHJvYnMsIDAuMiwgbXVzLCBzaWdtYXMsIGRmcykKCmRmICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT1vYnMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1zdGF0ZSkpCgpkcCA8LSBEaXJpY2hsZXRQcm9jZXNzR2F1c3NpYW4oc2NhbGUoZGYkb2JzKSwgZzBQcmlvcnMgPSBjKDAsIDEsIDIsIDEpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWxwaGFQcmlvcnMgPSBjKDAuMSwgNCkpCm5faXRlciA8LSA0MDAKZHBfMSA8LSBGaXQoZHAsIG5faXRlcikKZHBfMiA8LSBGaXQoZHAsIG5faXRlcikKCm5fY2x1c3RlcnNfMSA8LSBtYXBfZGJsKGRwXzEkbGFiZWxzQ2hhaW4sIH5tYXgoLikpCm5fY2x1c3RlcnNfMiA8LSBtYXBfZGJsKGRwXzIkbGFiZWxzQ2hhaW4sIH5tYXgoLikpCgp0aWJibGUobjE9bl9jbHVzdGVyc18xLCBuMj1uX2NsdXN0ZXJzXzIpICU+JSAKICBtdXRhdGUoaXRlcj1zZXFfYWxvbmcobjEpKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKC1pdGVyKSAlPiUgCiAgZ2dwbG90KGFlcyh4PWl0ZXIsIHk9dmFsdWUsIGNvbG91cj1uYW1lKSkgKwogIGdlb21fbGluZSgpCgpkcF8xJG51bWJlckNsdXN0ZXJzCmRwXzIkbnVtYmVyQ2x1c3RlcnMKZHBfMSRhbHBoYQoKZGYgJT4lIAogIG11dGF0ZShjbHVzdGVyPWRwXzIkY2x1c3RlckxhYmVscykgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPWNsdXN0ZXIpKQpgYGAKVHJ5IHRlbXBlcmluZyB0aGUgbGlrZWxpaG9vZApgYGB7cn0KTGlrZWxpaG9vZC50ZW1wZXJlZCA8LSBmdW5jdGlvbihtZG9iaiwgeCwgdGhldGEpewogIGJldGEgPC0gMQogIHZhbDEgPC0gMSAvIGJldGEgKiBkbm9ybSh4LCB0aGV0YVtbMV1dLCB0aGV0YVtbMl1dKV5iZXRhCiAgIyB2YWwyIDwtICgxIC8gKGJldGEgKyAxKSkgKiBpbnRlZ3JhdGUoZnVuY3Rpb24oeCkgZG5vcm0oeCwgdGhldGFbWzFdXSwgdGhldGFbWzJdXSleKGJldGEgKyAxKSwgLTIwLCAyMCkkdmFsdWUKICBhcy5udW1lcmljKHZhbDEpCn0KUHJpb3JEcmF3LnRlbXBlcmVkIDwtIGZ1bmN0aW9uKG1kb2JqLCBuPTEpewogIHRoZXRhIDwtIGxpc3QoKQogIHRoZXRhW1sxXV0gPSBhcnJheShybm9ybShuLCBtZG9iaiRwcmlvclBhcmFtZXRlcnNbMV0sIG1kb2JqJHByaW9yUGFyYW1ldGVyc1syXSksIGRpbT1jKDEsIDEsIG4pKQogIHRoZXRhW1syXV0gPSBhcnJheShyZXhwKG4sIG1kb2JqJHByaW9yUGFyYW1ldGVyc1szXSksIGRpbT1jKDEsIDEsIG4pKQogIHRoZXRhCn0KUHJpb3JEZW5zaXR5LnRlbXBlcmVkIDwtIGZ1bmN0aW9uKG1kb2JqLCB0aGV0YSl7CiAgcHJpb3JQYXJhbWV0ZXJzIDwtIG1kb2JqJHByaW9yUGFyYW1ldGVycwogIHRoZXRhRGVuc2l0eSA8LSBkbm9ybSh0aGV0YVtbMV1dLCBwcmlvclBhcmFtZXRlcnNbMV0sIHByaW9yUGFyYW1ldGVyc1syXSkKICB0aGV0YURlbnNpdHkgPC0gdGhldGFEZW5zaXR5ICogZGV4cCh0aGV0YVtbMl1dLCBwcmlvclBhcmFtZXRlcnNbM10pCiAgYXMubnVtZXJpYyh0aGV0YURlbnNpdHkpCn0KTWhQYXJhbWV0ZXJQcm9wb3NhbC50ZW1wZXJlZCA8LSBmdW5jdGlvbihtZG9iaiwgb2xkUGFyYW1zKXsKICBtaFN0ZXBTaXplIDwtIG1kb2JqJG1oU3RlcFNpemUKICBuZXdQYXJhbXMgPC0gb2xkUGFyYW1zCiAgbmV3UGFyYW1zW1sxXV0gPC0gb2xkUGFyYW1zW1sxXV0gKyBtaFN0ZXBTaXplWzFdKnJub3JtKDEpCiAgbmV3UGFyYW1zW1syXV0gPC0gYWJzKG9sZFBhcmFtc1tbMl1dICsgbWhTdGVwU2l6ZVsyXSpybm9ybSgxKSkKICBuZXdQYXJhbXMKfQp0ZW1wZXJlZE1kIDwtIE1peGluZ0Rpc3RyaWJ1dGlvbihkaXN0cmlidXRpb24gPSAidGVtcGVyZWQiLAogIHByaW9yUGFyYW1ldGVycyA9IGMoMCwgNSwgMC4wMSksCiAgY29uanVnYXRlID0gIm5vbmNvbmp1Z2F0ZSIsCiAgbWhTdGVwU2l6ZSA9IGMoMSwgMC4xKQopCmRwIDwtIERpcmljaGxldFByb2Nlc3NDcmVhdGUoc2NhbGUoZGYkb2JzKSwgdGVtcGVyZWRNZCkKZHAgPC0gSW5pdGlhbGlzZShkcCwgbnVtSW5pdGlhbENsdXN0ZXJzID0gMykKZHAgPC0gRml0KGRwLCAyMDApCgpkZiAlPiUgCiAgbXV0YXRlKGNsdXN0ZXI9ZHAkY2x1c3RlckxhYmVscykgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPWNsdXN0ZXIpKQpgYGAKCgoKVHJ5IGZpdHRpbmcgRFAgdG8gZmFrZSAyIG1lYW4gZGF0YQpgYGB7cn0KbiA8LSAyMDAKeCA8LSB2ZWN0b3IobGVuZ3RoID0gbikKZm9yKGkgaW4gc2VxX2Fsb25nKHgpKSB7CiAgdSA8LSBydW5pZigxKQogIGlmKHUgPiAwLjUpCiAgICB4W2ldIDwtIHJub3JtKDEsIDAsIDEpCiAgZWxzZQogICAgeFtpXSA8LSBybm9ybSgxLCAxMCwgMSkKfQp4IDwtIHNjYWxlKHgpCgpkcCA8LSBEaXJpY2hsZXRQcm9jZXNzR2F1c3NpYW4oeCwgZzBQcmlvcnMgPSBjKDAsIDEsIDIsIDEpKQpuX2l0ZXIgPC0gMTAwMApkcF8xIDwtIEZpdChkcCwgbl9pdGVyKQpkcF8yIDwtIEZpdChkcCwgbl9pdGVyKQoKbl9jbHVzdGVyc18xIDwtIG1hcF9kYmwoZHBfMSRsYWJlbHNDaGFpbiwgfm1heCguKSkKbl9jbHVzdGVyc18yIDwtIG1hcF9kYmwoZHBfMiRsYWJlbHNDaGFpbiwgfm1heCguKSkKCnRpYmJsZShuMT1uX2NsdXN0ZXJzXzEsIG4yPW5fY2x1c3RlcnNfMikgJT4lIAogIG11dGF0ZShpdGVyPXNlcV9hbG9uZyhuMSkpICU+JSAKICBwaXZvdF9sb25nZXIoLWl0ZXIpICU+JSAKICBnZ3Bsb3QoYWVzKHg9aXRlciwgeT12YWx1ZSwgY29sb3VyPW5hbWUpKSArCiAgZ2VvbV9saW5lKCkKCmhpc3Qobl9jbHVzdGVyc18yKQpgYGAKCgoKIyBSZWFsIGRhdGEKYGBge3J9CmRmIDwtIHJlYWQuY3N2KCIuLi8uLi8uLi8uLi9Eb3dubG9hZHMvZG9pXzEwLjUwNjFfZHJ5YWQuZHI3c3F2OXY5X192My9FbGVwaGFudHNEYXRhX2Fuby5jc3YiLHNlcD0iOyIpCgpkZiAlPiUgCiAgZmlsdGVyKGlkPT0xKSAlPiUgCiAgbXV0YXRlKHRpbWU9c2VxX2Fsb25nKGlkKSkgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PURpc3QpKSArCiAgZ2VvbV9saW5lKCkKYGBgCgpGaXQgRFAgdG8gZGF0YQpgYGB7cn0KbGlicmFyeShkaXJpY2hsZXRwcm9jZXNzKQoKZGZfMSA8LSBkZiAlPiUgCiAgZmlsdGVyKGlkPT0xKSAlPiUgCiAgc2xpY2UoMToyMDApCmRwIDwtIERpcmljaGxldFByb2Nlc3NHYXVzc2lhbihkZl8xJERpc3QsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhbHBoYVByaW9ycyA9IGMoMC4wMSwgNCkpCm5faXRlciA8LSA0MDAKZHBfMSA8LSBGaXQoZHAsIG5faXRlcikKZHBfMiA8LSBGaXQoZHAsIG5faXRlcikKCm5fY2x1c3RlcnNfMSA8LSBtYXBfZGJsKGRwXzEkbGFiZWxzQ2hhaW4sIH5tYXgoLikpCm5fY2x1c3RlcnNfMiA8LSBtYXBfZGJsKGRwXzIkbGFiZWxzQ2hhaW4sIH5tYXgoLikpCgp0aWJibGUobjE9bl9jbHVzdGVyc18xLCBuMj1uX2NsdXN0ZXJzXzIpICU+JSAKICBtdXRhdGUoaXRlcj1zZXFfYWxvbmcobjEpKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKC1pdGVyKSAlPiUgCiAgZ2dwbG90KGFlcyh4PWl0ZXIsIHk9dmFsdWUsIGNvbG91cj1uYW1lKSkgKwogIGdlb21fbGluZSgpCgpkcF8xJG51bWJlckNsdXN0ZXJzCmRwXzIkbnVtYmVyQ2x1c3RlcnMKZHBfMSRhbHBoYQpgYGAKClRyeSB3aXRoIG1vcmUgYW5nbGUgZGF0YSB0b28KYGBge3J9CmRmX3MgPC0gZGYgJT4lIAogIGZpbHRlcihpZD09MSkgJT4lIAogIHNsaWNlKDE6MjAwKSAlPiUgCiAgc2VsZWN0KERpc3QsIEFuZ2xlKSAlPiUgCiAgc2NhbGUoKQoKZGZfcyAlPiUgCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICBtdXRhdGUodGltZT1zZXFfYWxvbmcoRGlzdCkpICU+JSAKICBwaXZvdF9sb25nZXIoLXRpbWUpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT12YWx1ZSkpICsKICBnZW9tX2xpbmUoKSArCiAgZmFjZXRfd3JhcCh+bmFtZSkKCmRwXzEgPC0gRGlyaWNobGV0UHJvY2Vzc012bm9ybWFsKGRmX3MsIG51bUluaXRpYWxDbHVzdGVycyA9IDIsIGFscGhhUHJpb3JzID0gYygwLjAwMSwgNCkpCmRwXzIgPC0gRGlyaWNobGV0UHJvY2Vzc012bm9ybWFsKGRmX3MsIG51bUluaXRpYWxDbHVzdGVycyA9IDIsIGFscGhhUHJpb3JzID0gYygwLjAwMSwgNCkpCm5faXRlciA8LSAyMDAKZHBfMSA8LSBGaXQoZHBfMSwgbl9pdGVyKQpkcF8yIDwtIEZpdChkcF8yLCBuX2l0ZXIpCgpuX2NsdXN0ZXJzXzEgPC0gbWFwX2RibChkcF8xJGxhYmVsc0NoYWluLCB+bWF4KC4pKQpuX2NsdXN0ZXJzXzIgPC0gbWFwX2RibChkcF8yJGxhYmVsc0NoYWluLCB+bWF4KC4pKQoKdGliYmxlKG4xPW5fY2x1c3RlcnNfMSwgbjI9bl9jbHVzdGVyc18yKSAlPiUgCiAgbXV0YXRlKGl0ZXI9c2VxX2Fsb25nKG4xKSkgJT4lIAogIHBpdm90X2xvbmdlcigtaXRlcikgJT4lIAogIGdncGxvdChhZXMoeD1pdGVyLCB5PXZhbHVlLCBjb2xvdXI9bmFtZSkpICsKICBnZW9tX2xpbmUoKQpgYGAKTG9vayBhdCBjbHVzdGVycwpgYGB7cn0KcGxvdChkcF8xKQoKY2x1c3Rlcl9sYWJlbHMgPC0gZHBfMiRjbHVzdGVyTGFiZWxzCmRmX3MgJT4lIAogIGFzLmRhdGEuZnJhbWUoKSAlPiUgCiAgbXV0YXRlKGNsdXN0ZXI9Y2x1c3Rlcl9sYWJlbHMpCiAgcGl2b3RfbG9uZ2VyKC1pdGVyKSAlPiUgCiAgZ2dwbG90KGFlcyh4PWl0ZXIsIHk9dmFsdWUsIGNvbG91cj1uYW1lKSkgKwogIGdlb21fbGluZSgpCmBgYAoKIyBEUE0gb24gTWF0aGVtYXRpY2EgSE1NIGRhdGEKYGBge3J9CnggPC0gcmVhZC5jc3YoIi4uL2RhdGEvaG1tX3Rlc3QuY3N2IiwgaGVhZGVyID0gRkFMU0UpCgp4IDwtIHggJT4lIAogIHNsaWNlKDE6bnJvdyh4KSkKCmRwXzEgPC0gRGlyaWNobGV0UHJvY2Vzc0dhdXNzaWFuKAogIHNjYWxlKHgkVjIpLCBhbHBoYVByaW9ycyA9IGMoNCwgNCkpCmRwXzIgPC0gRGlyaWNobGV0UHJvY2Vzc0dhdXNzaWFuKAogIHNjYWxlKHgkVjIpLCBhbHBoYVByaW9ycyA9IGMoNCwgNCkpCmRwXzEgPC0gSW5pdGlhbGlzZShkcF8xLCBudW1Jbml0aWFsQ2x1c3RlcnMgPSAxKQpkcF8yIDwtIEluaXRpYWxpc2UoZHBfMiwgbnVtSW5pdGlhbENsdXN0ZXJzID0gMTApCgpuX2l0ZXIgPC0gMjAwMApkcF8xIDwtIEZpdChkcCwgbl9pdGVyKQpkcF8yIDwtIEZpdChkcCwgbl9pdGVyKQoKbl9jbHVzdGVyc18xIDwtIG1hcF9kYmwoZHBfMSRsYWJlbHNDaGFpbiwgfm1heCguKSkKbl9jbHVzdGVyc18yIDwtIG1hcF9kYmwoZHBfMiRsYWJlbHNDaGFpbiwgfm1heCguKSkKCnRpYmJsZShuMT1uX2NsdXN0ZXJzXzEsIG4yPW5fY2x1c3RlcnNfMikgJT4lIAogIG11dGF0ZShpdGVyPXNlcV9hbG9uZyhuMSkpICU+JSAKICBwaXZvdF9sb25nZXIoLWl0ZXIpICU+JSAKICBnZ3Bsb3QoYWVzKHg9aXRlciwgeT12YWx1ZSwgY29sb3VyPW5hbWUpKSArCiAgZ2VvbV9saW5lKCkKCmRwXzEkbnVtYmVyQ2x1c3RlcnMKZHBfMiRudW1iZXJDbHVzdGVycwoKaGlzdChuX2NsdXN0ZXJzXzEpCgojIGZpbmQgTUFQCm1fMSA8LSB3aGljaC5tYXgoZHBfMSRsaWtlbGlob29kQ2hhaW4pCm5fY2x1c3RlcnNfMVttXzFdCm1fMiA8LSB3aGljaC5tYXgoZHBfMiRsaWtlbGlob29kQ2hhaW4pCm5fY2x1c3RlcnNfMlttXzJdCmBgYAoKCgoKCg==